
# 灰色预测模型GM(1,1)实现示例
# 参考同目录下README.md的建模流程

import numpy as np
import matplotlib.pyplot as plt



# 导入原始序列
X0 = [86977, 99215, 109655, 120333, 135823, 159878, 182321, 209407, 246619, 300670]

# --- 可容覆盖区间验证（级比检验） ---
n = len(X0)
lambda_min = np.exp(-2/(n+1))
lambda_max = np.exp(2/(n+1))
lambda_list = [X0[i-1]/X0[i] for i in range(1, n)]
print("级比序列：", lambda_list)
print(f"可容覆盖区间: [{lambda_min:.3f}, {lambda_max:.3f}]")
c = 0  # 平移常数，默认为0
if all(lambda_min <= l <= lambda_max for l in lambda_list):
    print("原始数据满足级比检验")
else:
    print("原始数据不满足级比检验，需进行数据变换")
    # --- 数据平移变换 ---
    # 通过遍历c寻找使所有级比落入区间的c
    found = False
    for test_c in range(0, 100000, 1000):
        test_lambda = [(X0[i-1]+test_c)/(X0[i]+test_c) for i in range(1, n)]
        if all(lambda_min <= l <= lambda_max for l in test_lambda):
            c = test_c
            found = True
            print(f"找到合适的平移常数c={c}")
            break
    if found:
        X0 = [x+c for x in X0]
        print("平移后的数据：", X0)
    else:
        print("未找到合适的c，请人工调整或检查数据")


# 1. 构建累加生成序列
X1 = [X0[0]]
add = X0[0] + X0[1]
X1.append(add)
i = 2
while i < len(X0):
    add = add + X0[i]
    X1.append(add)
    i += 1
print("X1", X1)  # 累加序列

# 2. 计算邻均值序列
M = []
j = 1
while j < len(X1):
    num = (X1[j] + X1[j - 1]) / 2
    M.append(num)
    j = j + 1
print("M", M)  # 邻均值序列


# 3. 构建数据矩阵和向量
Y = []
x_i = 0
while x_i < len(X0) - 1:
    x_i += 1
    Y.append(X0[x_i])
Y = np.asmatrix(Y).T  # 向量Y
Y.reshape(-1, 1)
print("Y", Y)
B = []
b = 0
while b < len(M):
    B.append(-M[b])
    b += 1
print("B:", B)
B = np.asmatrix(B)
B.reshape(-1, 1)
B = B.T
B_const = np.ones((len(B), 1))  # 构建常数列
B = np.hstack((B, B_const))     # 数据矩阵B
print("B_const", B_const)
print("B", B)


# 4. 使用最小二乘法估计参数
beta = np.linalg.inv(B.T.dot(B)).dot(B.T).dot(Y)  # [a, b]^T = (B^T B)^{-1} B^T Y
a = beta[0]
b = beta[1]
const = b / a
print(beta)  # 参数a和b
print(type(beta))


# 5. 代入数据进行预测（GM(1,1)模型公式）
F = [X0[0]]
k = 1
while k < len(X0) + 10:
    F.append((X0[0] - const) * np.exp(-a * k) + const)
    k += 1
print("F", F)  # 累加预测序列


# 6. 还原，得到原始序列的预测值
x_hat = [np.asmatrix(X0[0])]
g = 1
while g < len(X0) + 10:
    print(g)
    x_hat.append(F[g] - F[g - 1])  # x^(0)(k+1) = x^(1)(k+1) - x^(1)(k)
    g += 1
X0 = np.array(X0)
x_hat = np.array(x_hat)
print(x_hat)


# 设置时间序列（原始与预测）
t1 = range(1999, 2009)      # 原始数据年份
t2 = range(1999, 2019)      # 预测年份

# --- 恢复数据的平移变换 ---
# 如果前面进行了数据平移（加了常数c），此处需将预测结果还原
if c is not 0:
    X0 = X0 - c
    x_hat = x_hat.reshape(-1) - c
    print(f"已将数据恢复为原始状态，c={c}")


# 结果可视化
plt.plot(t1, X0, color='r', linestyle="--", label='true')      # 原始数据
plt.plot(t2, x_hat.reshape(-1), color='b', linestyle="--", label="predict") # 预测数据
plt.legend(loc='upper right')
plt.xlabel('year')
plt.ylabel('Profit')
plt.title('Profit prediction for company by Grey model')
plt.show()

